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Abstract 



For the ordinary differential equation (ODE) x{t) — f{t,x), x{0) = Xq, t > 0, x ^ R'^, 
assume / to be at least continuous in t and locally Lipshitz in x, and if necessary, several times 
continuously differentiable in t and x. We associate a conditioning function E{t) with each 
solution x{t) which captures the accumulation of global error in a numerical approximation in 
the following sense: if x(t; h) is an approximation derived from a single step method of time 
step h and order r then \\x{t] h) - x{t)\\ < K{E{t) + e)h'^ for < i < T, any e > 0, sufficiently 
small and a constant K > Q. 

Using techniques from the stability theory of differential equations, this paper gives condi- 
tions on x{t) for E{t) to be upper bounded Hnearly or by a constant for < > 0. More concretely, 
these techniques give constant or linear bounds on E(t) when x{t) is a trajectory of a dynamical 
system which falls into a stable, hyperbolic fixed point; or into a stable, hyperbolic cycle; or 
into a normally hyperbolic and contracting manifold with quasipcriodic flow on the manifold. 

1 Introduction 

For the system of ordinary differential equations x{t) = f{t,x), t > 0, x € W^, the initial value 
problem x(to) = xq, Iq > 0, may not have a solution, and when the solution exists it may not be 
unique. However, if f{t, x) is continuous in t and not only continuous but also locally Lipshitz in x, 
there is a unique solution of the initial value problem for to < * < *o + £) e > 0, which we denote by 
x{t;to,xo)- We make these assumptions about f{t,x) and only consider solutions x{t;tQ,XQ) that 
can be continued till t = co. Usually to = 0, and we denote these solutions by x(t;xo) or simply 



Even at this level of generality, one might ask how accurately x(t;xo) can be approximated by 
a numerical method. A reasonable first guess is to look at the d x d matrix 



for < t < oo; since E'{t) gives the sensitivity of x{t;xo) with respect to xq, one might hope to 
relate it to the accumulation of global error. However, numerical methods introduce discretization 
error not just at t = but at every time step of the integration, and therefore, E'{t) proves to be 
an insufficient concept. 

*Research at MSRI is supported in part by NSF grant DMS-9701755 



x{t). 





dxo 



The following conditioning function E(t) associated with x{t;xQ), or more briefly x{t), does 
capture the accumulation of global error: 

E{t) = sup [' ^^v{s)ds , (1.1) 



v{s) 



dx{s) 



with the supremum taken over continuous functions v : [0,t] — > R'^ with ||v(s)|| < 1 for < s < t. 
Unlike E'{t), E{t) takes into account the sensitivity of x{t) with respect to all x{s), < s < t. In 
the definitions of E{t) and E'{t) above, we have assumed f{t,x) to be continuously differentiable 
with respect to t and x. All vector norms in this paper are Euclidean norms and all matrix norms 
are the corresponding induced norms. 

Let x{t; xq] h) denote the approximation to x{t; xq) computed by a single step method of order 
r. The numerical method gives x{t;xo;h) at all positive integer multiples of h. For kh < t < 
{k + l)h, /c = 0, 1, 2, . . . , we define x(r; xq; h) by following the solution exactly from the initial point 
t = kh and x = x{kh; xq; h) till t = t. Therefore, x(t; xq; h) can be discontinuous only at i = kh, 
k = 1,2, . . . The magnitude of this discontinuity 

\\x{{k + l)h; Xq; h) — x{{k + l)h; kh,x{kh; xq; ^))|| 

is the local discretization error at t = {k + l)h. We write this discretization error as Kf:^ih^~^^Vk+i, 
where Vk+i € R"^ with = 1 and i^fc+i is a non-negative real number. Thus the direction of 

the local discretization error at the ith step is Vi and its magnitude is Kih'^'^^. 

Let us assume that the magnitude of the local discretization error is bounded above by Kh^~^^ 
for a constant K > 0. This assumption can hold for example for a Runge-Kutta method of order r if 
f{t, x) is r + 1 times continuously differentiable with respect to t and x; we discuss this assumption 
further in Section 2. Then by Theorem given e > and T > 0, 

||x(t;xo) -x(t;xo;/i)|| < {E{t) + e)Kh'' (1.2) 

for < t < T and sufficiently small h. In this bound on the global error, E{t), which is given by 
(|l . l| ) , is independent of details of the numerical method. Further, a bound like (|l.2D will not hold 
for any real-valued function of t which is strictly less than E{t) at some value of t. It is because of 
these reasons that we call E{t) a conditioning function. 

Let us now draw an analogy between E{t) and the absolute condition number of a multivariate 
function ^(x). The formula (see p^) 

\\g{x + 6x) - g{x)\\ 
hm sup 

5^0+ Sx<5 

is the analogue of ( |l.lD ; this absolute condition number measures the sensitivity of g{x) with respect 
to small changes in x. For a stable numerical method, the error in evaluating g{x) is governed by 
this absolute condition number in a manner similar to the dependence of global error on E{t) given 

by (0). 

Sections 2,3, and 4 define E{t) and develop its properties. We do not begin with (1.1) as the 



definition of E{t). Section 2 introduces a model of discretization errors which is similar to but more 
general than the model of Stuart and Humphries |2H|. Section 3 defines E{t) using this model. The 



expression for E{t) in (1.1) is derived in Corollary 6.3 of Section 6. Let us mention the formal 



similarity of ( |1.1D to bounds on the global error derived by Iserles and Soderlind [16| using the 
Alexseev-Grobner lemma. The results in this half of the paper argue that E{t) is the appropriate 
vehicle for a study of global errors. 
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From Section 5 onwards, we undertake the task of relating E{t) to stability properties of x{t; xq). 
The relation of E{t) to the accumulation of global error as t increases is clear enough from ( |l.2D . 
If E{t) is bounded by a constant or linearly or by a polynomial of low degree in t, the accurate 
approximation of the trajectory x{t; xq) can be considered a tractable problem; but if E{t) increases 
exponentially in t accurate approximation of x{t; xq) is pretty intractable. The standard technique 
of bounding global errors using the Lipshitz constant is used in convergence proofs of numerical 
methods @. But the bounds on E{t) obtained this way increase exponentially in t and are of 
hardly any other use. 

We upper bound E{t) by making stability assumptions on x{t;xo). Sections 5 and 6 derive 
linear and constant bounds on E{t) by making stability assumptions on x{t; xq). This connection to 
stability is somewhat subtle; there are exponentially stable examples with exponentially increasing 
E(t). However, the work of researchers who followed Lyapunov allows us to clarify and circumvent 
the difficulties in relating E{t) to stability theory. See Table |l| in Section 6 for a summary. 

Let us now make some prefatory remarks about the stability theory of ordinary differential 
equations. The theory of stability of ordinary differential equations was initiated by A.M. Lyapunov 
in 1892 One stream of research which emanates from this remarkable work is about first 

approximations. Let x{t) = x{t;xo) be a solution of = f{t,x). If the perturbation y{t) is such 
that x{t) + y{t) is also a solution of the same equation, then y(t) satisfies 

= fit, y + xit)) - fit, xit)) = Fit, y), 

where clearly Fit, 0) = 0. The solution is stable in the sense of Lyapunov if for every e > 
there exists a. 5 > Q such that ||y(0)|| < 5 implies ||y(t)|| < e for t > 0. Now (see Q), xit) is 
stable if and only if the zero solution of = y) is stable. Thus it is enough to look at zero 
solutions of systems of the form = y), with -F(t, 0) = 0. 

Lyapunov pointed out an other possible simplification. Since Fit, 0) = 0, if Fit, y) is assumed 
to be continuously differentiable in y, -F(t, y) = Ait)y + o(||y||) as y ^ 0. Here, 

dFit,y) 



Ait) 



dy 
dfit,x) 



(1-3) 

ox x=x(t;xo) 

Lyapunov argued that perhaps stability properties of the zero solution of the linear first approxi- 
mation y(t) = Ait)y might imply stability properties of the zero solution of the nonlinear system 
y(t) = F(t, y). This program has been carried out by Lyapunov, E. Cotton (1911), O. Perron 
(1928), I.G. Petrovskii (1934), R. Bellman (1953), and others. 

The other stream of research, which also originated with Lyapunov, uses Lyapunov functions 
Vit, y). If every solution y(t; yo) of y(t) = F(t, y), i > 0, had the property that ||y(t; yo)|| decreases 
as t increases, stability of the zero solution would be easy to infer. However, this property does 
not hold even for stable, linear systems of the form y(i) = Ay, where ^ is a constant matrix. 
Lyapunov's idea was to use a non-negative, real- valued function Vit,y), instead of the norm, such 
that y(t, y(i)) decreases as t increases. This F(t, y) is required to be related to ||y|| in a uniform 
way for i > 0; precise details depend upon the concept of stability. After Lyapunov, this line 
was greatly developed by Soviet researchers, including Persidskii, Malkin, Krasovskii, and others, 
beginning in the 1920s. It was taken up by researchers in the west, including J.L. Massera, T. 
Yoshizawa, and others, in the 1950s. 

Parts of Sansone and Conti |22] and Hale |1C] are excellent introductions to the theory of 
stability of ODEs. There are numerous advanced works; of them we refer to Bellman Malkin 
|T8|, and Yoshizawa [p^l. 
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Obtaining detailed stability information about solutions of ODEs can be far from trivial. For 
this reason, applying the theory in Sections 5 and 6 to concrete examples is not a simple matter. 
In Section 7, we give three applications to dynamical systems. We derive constant or linear upper 
bounds on E{t) for trajectories falling into stable, hyperbolic fixed points, or into stable, hyperbolic 
cycles, or into a normally hyperbolic and contracting manifold, with the flow on the manifold being 
quasiper iodic. These three applications involve the Hartman-Grobman theorem, convergence in 
phase results for stable cycles, and results of Pugh, Shub, and others about normally hyperbolic 
flows, respectively. The first of these applications was covered by Stuart and Humphries p5[| . The 
second application is a significant extension of a result due to Cano and Sanz-Serna Q. The third 
application appears to be entirely new. 

Let us mention the similarity of our analysis to the asymptotic analysis of global error of 
Henrici [12| |13] and several following researchers, including Gragg j^. Theorem |6.2| , in particular. 



is implicit in Henrici's work. Beginning with Dahlquist [^, it has been known that the use of 
one-sided Lipshitz constants can sometimes give meaningful bounds on the global error. We show 
in Section 8 that the use of one-sided Lipshitz constants fits naturally into our framework. Section 
8 also considers variable time stepping, multistep methods, and other issues. 

To summarize briefly, this paper consists of three parts. The first part, Sections 2, 3, and 4, 
derives a conditioning function E{t) and associates it with the global error in numerically approxi- 
mating the solution of an ordinary differential equation. The second part. Sections 5 and 6, relates 
E{t) to the stability theory of ordinary differential equations. The third part, Section 7, applies 
this theory to dynamical systems. 

2 A Model for Discretization Errors 

The model of discretization error which we now present is close to the way discretization errors 
are made by single step methods with constant step sizes. Stuart and Humphries [p5|| model 
discretization errors of single step methods in a similar manner. 

Let a{h) be a continuous, strictly increasing function of h for h > 0. Assume also that a(0) = 0. 
Then an approximation Xa{t; xq; h) to x{t; xq) is defined as follows: 

Xq(0;xo; h) = xq 

Xainh; Xq] h) = x{nh; {n — l)h, Xa{{n — l)h; Xq] h)) + ha{h)vn n>l, (2.1) 

where Vn G R"^ can be any vector with ||w„|| < 1. In words, the approximate solution at t = nh, 
n > 1, is obtained by exactly propagating the point Xa((n — l)h;xo;h) at i = (n — l)h under 
x{t) = f{t, x) till t = nh, and then adding the discretization error or discontinuity ha{h)vn, where 
ll^nll ^ 1- For (n — l)h < t < nh, 

Xa{t; Xq; h) = x{t; {n — l)h, Xa{{n — l)h; xq; h)). (2-2) 

Since Vn can be any vector with \\vn\\ < 1, this actually defines a whole family of approximate 
solutions which we denote by Xa{xQ; h). The exact solution x{t) is the only member of this family 
which is continuous. Figure |l| gives an example of an approximate solution. 

We now address how single step numerical methods are related to approximations from the 
family X{xQ;h). In our description of single step methods, the discretization error at every step 
was taken to be Kih^^^Vi. We now assume that the Ki are bounded by a constant K which does 



not depend upon h or i. This can be proven in some circumstances; see [25|. Besides, if there 



is no such K, the numerical method will not in practice behave as if it were of order r. Thus 
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Figure 1: The thick hnes show an exact solution of the equation x{t) = sint + (cost)x and an 
approximation to it with h = 0.5 and a{h) = 2h for. 



when the order of accuracy of the numerical method for approximating x{t;xo) is r, we can take 
a{h) = Kh^, and there will always be an approximation in the family Xa{xo; h) which is the same 
as the trajectory of the numerical method. 

The notations we have established so far will be adhered to in the rest of the paper. For the 
system, x{t) = f{t,x), f{t,0) is not necessarily zero. The solution with x(to) = xq is denoted by 
x{t;to,xo), and by x{t;xo) when to = 0. The approximations to x{t;xo) obtained as in (2.1) and 
(|2.2| ) are denoted by Xa{t;xo;h). The family of approximations is X{xQ;h). The same notations 
apply for y{t) = F{t,y) with the xs changed to ys. But here F{t,0) = 0, and ya{t;h), t > 0, is 
an approximation to the zero solution, and Ya{h) is a collection of those approximations. When 
we speak of a linear first approximation, it is always obtained as in ( |1.3D . Sometimes we omit the 
word linear. This same equation is sometimes called the equation of first variation or simply the 
linearization. Let us note that when we speak of stability, it is always stability in the sense of 
Lyapunov and his followers; in particular, it is not numerical stability in the sense of Dahlquist. 



3 Definition of Ea{t) 

Let us recall that a{h) is assumed to be a strictly increasing, continuous function of h for h > 
which is zero at zero; for example, a{h) can be Kh^ for a positive integer r. Assume the functions 
f{t,x) and F{t,x) to be defined for < t < oo and x € R"^, to be continuous in t, and locally 
Lipshitz in x. Besides, F{t,0) = 0. Assume the solution x{t;xQ) of the initial value problem 
x{t) = f{t,x), x(0) = xo, to be continuable till t = oo. Only some of these assumptions are 
restated in theorems that follow. 

The global error ea{t;xQ; h) is defined as follows: 

ea{t;xo;h)= sup \\xa{t;xo;h) - x{t;xo)\\. 
ex 

We later use this to define Ea{t). 
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Instead of saying that Xa{t; xq; h) exists for < t < T, we say that Xa{t; xq; h) can be continued 
till t = T. We say that every approximation Xa{t; xq; h) can be continued till T if every Xa{t; xq; h) 
exists for < t < T with any allowed choice of discontinuities at t = kh. Lemma 3.1 introduces 
ho(T,r) and L(T,r). 

Lemma 3.1. Assume, as usual, that f{t,x) is continuous in t, t > 0, and locally Lipshitz in x, 
X S . Let x{t;xQ), t > 0, be the unique solution of the initial value problem x{t) = f{t,x), 
x(0) = Xq. Then, 

(i) There exists a constant L{T, r) > such that 

\\f{t,xi)- f{t,X2)\\ < L(t,r)||xi -X2II 
for < t <T and \\xi — x{t; xq) \\ < r, for any T > 0, r > 0, and i = 1,2. 

(ii) There exists a constant hQ{T, r) such that for < h < hQ{T, r) every approximation Xa(t; xq; h) 

can be continued till t = T, and satisfies \\xa{t;xQ;h) — 2;(t;xo|| < r, for < t < T . 

Further, forO<h< ho{T,r), ea{t;xo;h) < te^('^''^)*a(/i). 



Proof. Similar to proof of Theorem 3.4.6 in Stuart and Humphries |25|. Similar estimates using 
the Lipshitz constant are found in Q and other places. □ 

For the zero solution of y{t) = F{t, y), y(0) = 0, where F{t, 0) = 0, ea{t; h) is defined as follows: 

ea{t]h)= sup ||ya,(t;/i)||, (3.1) 

where Ya{h) is the family of approximations to the zero solution with time step h. 

Proposition 3.2. Assume h < /io(r, r) for r > 0, T > 0, and that < t < T . The ea{t;xQ;h) 
of the solution x{t;xo) of x{t) = f{t,x), x(0) = xq, and the ea{t;h) of the zero solution of y{t) = 
F{t, y), y(0) = 0, where F{t, y) = f{t, y + x{t; xq)) — f{t, x{t; xq)), are the same. 

Proof. It is enough to show that members Xa(t; xq; h) of Xa{xo; h) and members ya{t] h) of Ya{h) 
can be matched so that Xait; xq; h) = x(t; xq) + ya(t; h). 
It is well known (see Q) that 

x{t + r; t; x{t; xq) + S) = x{t + r; xq) + y{t + r; t, 5) 

for T > 0. Hence, if Xa{kh; xq; h) = x{kh; xq) + ya{kh; h) then Xa{t; xq; h) = x{t; xq) + ijait', h) for 
kh < t < {k + l)h. Further, the discontinuities of Xa{t; xq; h) and yaif-'-, h) at t = kh, where A: is a 
positive integer, can be exactly matched. Therefore, the approximate solutions can be matched as 
desired. □ 

The Ea{t) which corresponds to the zero solution of y{t) = F{t,y), y(0) = 0, i > 0, is defined 
as follows: 

i?„(t) = limsup^4#- (3.2) 

By Lemma |T], Ea{t) < te^\ where L = L{T, r) for some T > t and r > 0. For a nonzero solution 
x{t; Xq) of x{t) = f{t, x), x{Q) = Xq, t > 0, Ea{t) is defined to be the same as the Ea{t) for the zero 
solution of y{t) = f{t, y + x{t; xq)) - f{t, x{t; xq)), y(0) = 0, i > 0. 
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As in the stability theory of ODEs, we can and do confine ourselves to an analysis of the zero 
solution of y{t) = F{t,y) without any loss of generality. Prom here on, assume L(T,r) to be such 
that 

\\F{t,yi)-F{t,y2)\\<L{t,r)\\yi-y2\\ 

for < t < T, \\yi\\ < r, \\y2\\ < r, where T > and r > 0. Also, assume hQ{T,r) to be such 
that for < h < hQ{T,r) every approximation ya{t;h) can be continued till t = T and satisfies 
\\ya{t;h)\\ < r. 



4 Properties of Ea{t) 



This section is devoted to properties of Ea{t). We show that Ea(t) is continuous (Theorem 4.3 ) 



and independent of the choice of a{h) (Theorem 4.4) so that it can be written as E{t). Theorem 



4.5 and its corollary relate E(t) to the accumulation of global error when x{t; xq) is approximated. 

In the proofs in this section, we repeatedly use Theorem 10.1 of Q. That theorem allows us to 
bound the divergence of two solutions of a differential equation using a Lipshitz constant. 

A sequence of inequalities are often combined in a way we now describe. Let eg = 0, and 
Cj < fci-i + rj_i for 1 < i < n. Then, we have, 

e„< (r-Vo + r-Vi + ---+r„_i), 

assuming / > 0. In fact, most of the time rj will all be equal. In that situation, e„ < ro(l + / + 
• • • + /"~^). And when / > 1, we can get e„ < f^^^{rQ + • • • + r„_i), or e„, < n/"~^ro if the are 
all equal. 

Lemma 4.1. Let h < h(){t + s, r), r > 0. Then 

ea{t + s;h) < e^^ea{t; h) + {s + h)e^^a{h), 
where L = L{t + s, r) and s > 0. 

Proof. Let ya{t;h) be an approximation to the zero solution of y{t) = F{t,y), 7/(0) = 0. 

Let t + hi be the first multiple of h after t, t + s — hr the last multiple before t + s, and let 
t + s — hr = t + hi + nh. Using Theorem 10.1 of and taking into account the discontinuities at 
multiples of h, we have 

\\ya{t + hi;h)\\ < e^^'\\ya{t;h)\\ + ha{h), 
\\ya{t + hi + kh; h)\\ < e'^^||yQ,(t + hi + {k - l)h; h)\\ + ha{h), 1 < k < n, 
\\ya{t + s;h)\\ < e^''^\\ya{t + nh;h)\\. 

Combining these inequalities, we get 

\\ya{t + s;h)\\ < e^'\\ya{t-h)\\ + {s + h)e^'a{h). 



The proof is easy to complete using (|3.l|) . □ 



Lemma 4.2. With the same assumptions about h and L as in the previous lemma, 

ea{t + s;h) > e'^^eait; h), 

where s > 0. 
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Proof. Let ya{T;h), < r < t, be an approximate solution. Continue it from T = t to T = t + s 
exactly as a solution of y{t) = F{t,y). Then, 

\\ya{t + s;h)\\ > e-^'\\ya{t;h)\\. 

The proof is now easy to complete using ( p.lD . □ 

Theorem 4.3. The Ea{t) of the zero solution of y{t) = F(t,y), y(0) = 0, is continuous for t >0. 



Proof. From Lemmas 4.1 and 4.2, we get 

ea{t - 6; h)e-^^ < e„(t; h) < e^^ea{t -6;h)+ e^\5 + h)a{h), 
when < t — (5 < and when t <t + 5, 

{ea{t + 6;h)- e^\5 + h)a{h))e-^^ < ea{t; h) < e^^ea{t + d;h). 
Divide these two inequalities by a{h) and use (|3.2|) to deduce continuity of Ea(t). 



□ 



Theorem 4.4. The Ea{t) of the zero solution of y{t) = F{t,y), y{0) = 0, is the same for any 
a{h), with a(0) = and a{h) continuous and strictly increasing for h >0. 

Proof. Let P{h) be another function like a{h), which is continuously increasing for h >0 and zero 
at zero. We will show that Ep[t) = Ea{t). 

Take L = L{t, r) for some r > 0. Take ha and hp small enough that the approximate solutions 
in the families Ya{h) and Yi3{h) are not only continuable till t but stay within a radius r of zero. 
For e small enough, there exist unique ha and hp such that a{ha) = Pihp) = e. Assume this 
relationship between ha, hp, and e. Then /i^ — > and /i^ — > as e ^ 0. 

It is enough to show that 

\ea{t;ha) - ep{t;hp)\ 
lim = U. 

In fact, it is enough to show that for every approximate solution ya{t', ha) there is an approximate 
solution yp{t;hp) such that 

\\ya{t; ha) - yp{t; hp)\\ < ef{t)g{e), (4.1) 

where f{t) is finite and (7(e) — > as e ^ 0. It will also have to be shown that a given yp{t; hp) 
can be approximated by some ya{t;ha) in the same manner; but the proof of this is gotten by 
transposing a and p. 

Given a ya{t; ha), we construct a yp{t; hp) as follows. yp{T; hp) is determined by the differential 
equation y{t) = F{t,y), except when r = php, p being a positive integer. At points r = php, 
introduce a discontinuity of magnitude smaller than or equal to hp[5{hp) in such a way as to get as 
close to ya{php; ha) as possible. 

Let hp = {m + f)ha, m > being an integer, and < / < 1. Assume that the number of 
integer multiples of ha in {php, (p + l)hp] is kp. Let Cp = \\ya{php; ha) - yp{php; hp)\\. 
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Then cq = and 



ep+i < e'^^^Cp + kpho,a{h^)e^f<^ - hpP{hp) 
= e'^^^Cp + kphaee^^^ - h^e. 

The first two terms can be derived from Theorem 10.1 of |^]. The third term is because of the 
discontinuity introduced into ^/^(r; hp) at r = phf^. If this bound on e^+i is negative, we can take 
ep_|_i = 0. If n = L^J) s-iid there are kn multiples of ha in then 



\ya{t; ha) - y/sit; hp)\\ < e^'^^'en + knhaa{ha) 



e 



where hr = t — nhp. 

Let p be the highest among {0, 1, . . . ,n} such that Cp = 0. Combining the inequalities for e,, 

p < i < n, we get 

en < e^("-P)'^'3(/i„ee''''^(fcp + • • • + kn-i) - {n-p)hpe). 

But kp + ■ ■ ■ + kn-i being the number of multiples of ha in [php^nhp] is bounded above by (n — 
p)hp/ha + 1. Let t' = t — hy. Then, 

en < e^^'{{n - p)hpe^i'^e + hace^i'^ - (n - p)h(se) 
< max(t'e^*',e^*')((e^'3^ - I) + hae''''^)e. 

Finally, 

WVait; ha) - yp{t- hp)\\ < max(e^*, te^*)((e^'3^ - 1) + hae^^^ + A:„/i«)e. 

Since kn is the number of multiples of ha in (t — hr,t], knha < hp + ha- Thus the bound above 
satisfies all conditions required of and the proof is complete. □ 

From here on, we drop the subscript a from Ea{t), ya{t'i h), Ya{h), Xa{t',xo, h), Xa{xo; h), and 
ea{t;h). 

Theorem 4.5. Let E{t) be associated with the zero solution of y{t) = F{t,y), y(0) = 0. Given 
T > and e > 0, there exists /iq > such that h < ho implies 

sup|iy(t; h)\\ = e{t; h) < {E{t) + e)a{h) 

forO<t<T. 

Further, if 9 < E[t) there exists an e > such that e{t; h) > [6 + e)a{h) for arbitrarily small h. 



Proof. Let e(t; h) = e(t; h)/a{h) - E{t). 

Take h < /iq, where /iq to begin with is smaller than ho{T,r) for some r > 0, and L = L{T,r). 
By definition of E{t) Q, 

limsupe(t,/i) = (4.2) 

for < t < T. 
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To prove the first part, given e > we find an /iq > such that h < implies e(t;h) < e for 
< t < T. In fact, it is enough to show that there exists an ht > 0, which depends on t, and an 
open neighbourhood of t in [0, T] such that for r in that neighbourhood and h < ht, e(r; h) < e. If 
the neighbourhoods of ti, . . . ,t„ give a finite cover of the compact interval [0, T], we can take ho 
to be the minimum of /i^ ^ , . . . , ht^ . 

Clearly, 

e{t + 6t; h) = e{t; h) + + ^) - M _ ^^^^ ^ _ ^^^^^ 

By (U), e(t;/i) < ei for /i smah enough. By continuity of E{t), \E{t + 5t) - E{t)\ < for \bt\ 
small enough. 

Using Lemmas [4.1| and E?^, and the bound iE^^a{K) on e(t; /i) from Lemma 3.1, we get 



a{h) 

for 5t > 0, and for 6t < 0, 

\e{t + 6t;h)-e{t;h)\ f rt,_Lst ,rtn „L5u 



a{h) 



< max(^te^*(e-^^* - l),te^*(l - e^^^) + (5t + /i)). 



Therefore, we can choose \6t\ and h small enough to ensure \e(t + 5;h) — e{t; h)\/a{h) < £3. We 
need only take ei + e2 + £3 < e to find the desired ht and the neighbourhood of t. 

The second half of the theorem is immediate from the definition of E(t). □ 

Corollary 4.6. Let E[t) he associated with the solution x(t;xo) of x{t) = f{t,x), x{0) = xq, 
where f{t,x) is continuous in t and locally Lipshitz in x in an open neighbourhood of the solution 
{t,x{t;xQ)), t >0. Given T > and e > 0, there exists an ho > such that 

\\x{t; xq] h) - x{t; xo)\\ < {E{t) + e)a{h) 

for < t < T , h < ho, and any approximation x{t; xq; h) in the family X{xq; h). 

So far, we have defined the notion of an approximate solution (see Figure ||), and an E{t) for 
every solution of an ODE. Theorem and its Corollary 4.6 relate E{t) to the accumulation of 
global error. In the upper bound {E[t) + e)a{h) for the global error, the details of the numerical 
method have been pushed into a{h); the conditioning function E{t) is something intrinsic to the 
solution x{t]XQ). Figure § shows E{t) for some simple examples. 

As shown in Figure |2|a, the bound on global error given by E{t) may be pessimistic in practice. 
Our model of approximations allows arbitrary discontinuities whose norms are bounded by ha{h), 
while for any given numerical method the discretization errors are fixed. Our analysis does not 
say how small an h is good enough for the bound on the global error to be governed by E{t) as in 



Theorem 4^. Will an /i in a practical computation be small enough? We answer this in Henrici's 
words [13|: usually if a stepsize is small enough to yield an accurate solution, it is also small 
enough that an asymptotic formula gives a correct indication of the size of the error. Numerical 
experiments in [0] and 0] lend ample support to this statement. 
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Figure 2: (a) The thick line is the E{t) for y{t) = -{2{t + 1) + {t + l)~^)y, y(0) = 1. The two 
thin hnes below it are global errors of forward and backward Euler divided by 8h; here h = .01. 
(b) E{t) for y(t) = Xy. The exact E{t) are obtained using Theorem 5.1 . 



5 E{t) for Linear Systems 

Our investigation of the relationship between E{t) and stability properties of the exact solution 
begins with the linear system y{t) = A{t)y, y(0) = 0. The relationship is not as simple as one 
might wish. There are both asymptotically stable examples with exponentially increasing E{t) and 
unstable examples with linearly bounded E{t). However, Theorems 5.5 and 5.6 give conditions for 
E{t) to be bounded by a constant or to be linearly bounded. Linear systems are an important 
class of problems by themselves. What is more, they can be used for understanding the E{t) of 
nonlinear systems. We will see this in Theorem |6.2| and in Section 6.2. 

In y{t) = A(t)y, the d x d matrix A{t) is assumed to be continuous for t > 0. For such linear 
systems, it is easy to show that y{t) is a linear function of y(0) for t > |22]. Thus y{t) = Y{t)y{0) 
for Y{t) G i?'^''^, where Y{t) is called the principal fundamental matrix of the linear system. The 
matrix Y{t) itself is continuous in t and always nonsingular. Moreover, y{t) = Y {t)Y~^ {s)y{s) for 
t > s. Although it does not appear to be directly related, let us mention the beautiful derivation 
of the Magnus series for Y{t) by Iserles and N0rsett |15]. 

For scalar linear systems y{t) = a{t)y, a{t) € -R, we have the following theorem. 

Theorem 5.1. The E{t) for the zero solution of y{t) = a{t)y, y(0) = 0, is given by 



where 



Jo 



g{t) = / a{T)dT. 
Jo 



Proof. The fundamental matrix, which is scalar in this situation, is given by Y(t) = e^^^\ Since 
Y{t) is always positive, the optimal choice of v{s) in Theorem 5.2 is v{s) = 1. □ 
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Let us now consider the concepts of stability put forward by Lyapunov [|^. The definitions 
will be stated for nonzero solutions x{t;xQ). 



Definition 1. The solution x{t;xo) of x{t) = f{t,x), x(0) = xq, is stable if given any e > there 
exists a 6 > such that ||xq — xo|| < S implies ||x(t;xQ) — x(t;xo)|j < e for t > 0. In fact, stability 
implies that given e > 0, there exists a 6{t) > for every r > such that ||x(r; Xq)—x{t; xo)\\ < 5{t) 
implies that ||x(t;xo) — x{t;xo)\\ < e for t > r. Let us emphasize that 6{t) can depend on r. 

Definition 2. The solution x{t;xQ) is asymptotically stable if given e > there exists a 5{t) > 
for every r > such that \\x'^ — x{t;xq)\\ < 6{t) implies not only that \\x{t;T,x'^) — x{t;xQ)\\ < e 
foi t > T but also that \\x{t;T,x'^) — x{t;xQ)\\ — > as t ^ oo. 

Implicit in the definitions is an assumption about the existence of solutions which begin near the 
solution x{t;xo). Obviously, asymptotic stability implies stability. For the scalar, linear problem 
y{t) = a{t)y, y(0) = yo, a necessary and sufficient condition for asymptotic stability is g{t) — > — oo 
as t ^ oo, where g(t) = Jq a{s)ds. However, the following examples show that both these concepts 
of stability are insufficient for bounding E{t). 

Example 5.1. We will show using the scalar, linear problem y(t) = a{t)y that even when a solution 
is asymptotically stable the E{t) associated with it can increase at an arbitrarily high rate. Given 
a rate r{t), consider a continuously differentiable function g(t), t > 0, g{0) = 0, such that 

1. g{t) < -t for all t > 0, 

2. ef e-9('Us > r{k) for A: = 1, 2, 3, . . . 

For the linear system, take a{t) = g'{t). The first condition ensures asymptotic stability of the zero 
solution, and the second condition implies E{k) > r{k) for positive integers k. Such a g{t) is easy 
to construct. Take g{k) = —k for k = 0, 1, 2, ... . For k — l<t<k^k>l, define g{t) so that 
g{t) < —t and 



This can be carried out for any continuous r{t), for example r{t) = e*. 

Example 5.2. On the other hand, there are unstable solutions with linearly bounded E{t). Con- 
sider the scalar, linear system y{t) = j^y, t >0. For this ODE, y{t) = {l + t)'^y{0) implying insta- 
bility of the zero solution for a > 0. Yet, for < a < 1 E{t), which is (l-a)-^(l+t)"((l+t)^-"-l), 
is linearly bounded. For a = 1, E{t) is (1 + t) log(l + t). 

Example 5.3. In Example |5.1| , \a{t)\ will be unbounded. Does asymptotic stability of = a{t)y, 
t > 0, imply a linear bound for E{t) if \a{t)\ is bounded? The answer is no; E{t) can still increase 
exponentially in t. We sketch the construction of a g{t) to show this. First take gi{t) = —t and 
g2{t) = — 2t. Take g(t) = g2{t) for < i < ti, and let g{t) increase monotonically till g{t2) = 51(^2) 
for t2 > ti, and then let g{t) decrease monotonically till gits) = 52(^3) for ^3 > t2. Repeat the same 
construction from t^ onwards with t4, ts, and t^ in place of ti, t2, ^3, and so on. The construction 
may be arranged so that 

1. if g{T) = 5-1 (r) then g{t) = g2{t) for fiT <t< f2T for any fixed < /i < /2 < 1, 




■k 
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2. \a{t)\ = \g'{t)\ is bounded. 

It is easy to check that E{t) > e~'^^'^{e^^^'^ — e^-^^"^), for r such that ^(r) = gi{T). Therefore, for 
/2 > 1/2, E{t) increases exponentially. Let us note that the linear system in this example has a 
negative Lyapunov exponent of —1. 

Notions of stability needed for bounding E{t) either by a constant or linearly are introduced 
after Theorem 5.2. Theorem 5.2 is comparable to Theorem 3.1 of Cano and Sanz-Serna [^]. 

Theorem 5.2. The E{t) of the zero solution of y{t) = A{t)y, y(0) = 0, is given by 



E{t) =sup j Y{t)Y~^^ 

v(s) Jo 



s)v{s)ds 



where the supremum is over all continuous functions v{s) : [0, t] —>■ R'^ with ||f (s) || < 1 for < s <t. 
As before, Y{t) is the principal fundamental matrix of y{t) = A{t)y and A{t) is continuous. 

Proof. As usual, let the discontinuity of y{t; h)att = kh be ha{h)vk. Let n = and hr = t — nh. 

Then 

y{kh; h) = Y{kh)Y-\{k - l)h; h)y{{k - l)h; h) + ha{h)vk, 
for 1 < A; < n, and y{t; h) = Y{t)Y^'^{nh)y{nh; h). Combine these equalities to get. 



y{t-h) = ^Y{t)Y-\kh)haih)vk, 



(5.1) 



k=l 



This expression is also used in the proof of Theorem 6.2 
To prove, 

ft 



sup 

v{s) 



[ Y{t)Y~h 
Jo 



s)v{s)ds < E{t), 



it is enough to find a y(t; h) for every v{s) such that 



'Y{t)Y-\sHs)ds=y^ + v{h), 



(5.2) 



(5.3) 





where r/(/i) — > as /i — > 0. Take the discontinuity of y{t; h) at t = kh to be ha{h)v{kh). Use i f>.l\ ) 
and it is apparent that y(t; h)/a{h) approximates the Riemann integral fQY{t)Y~^{s)ds. The 
standard convergence theorem for Riemann integrals gives r](h) ^ as /i ^ 0. 

To show ( ^ ) in the reverse direction, it is enough to find a v{s) for a given y{t; h) so that ( ^.31 ) 
holds. First, consider a discontinuous v{s) defined by v{s) = for kh < s < {k + l)h. Clearly, 
\\v{s)\\ < 1. Then, 







'Y{t)Y-\s)i}{s)ds = + 6{h)t, 

a[h) 



(5.4) 



where 6(h) = supQ^g^f\\Y{t)Y~^{s + h) — Y{t)Y~'^{s)\\. Since Y{s) is continuous, S{h) ^ as 
/i — >• 0. Lusin's theorem I^T]] is a basic result for approximating measurable functions by continuous 
functions. It guarantees a continuous f(s) such that ||u(s)|| < 11^(5)11 and //(i'(s) / v{s)) < h, 
where ^ is the Lebesgue measure. If M = supo<s<( \\Y{t)Y-\s)\\, then 

Y{t)Y-^{s)v{s)ds = / Y{t)Y-'^{s)d{s)ds + 2mh, \m\<2M. (5.5) 
Jo 

Together, (5.4) and ( |5.5| ) complete the proof. □ 
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Corollary 5.3. 



Eit) < [ \\Y{t)Y~\s)\\ds. 
Jo 



Corollary 5.4. For < 6 < t , 

[■5 

E{t) > II / Y{t)Y~\s)ds\ 
Jo 



The definitions of uniform stability and uniform asymptotic stability that follow seem to have 
been introduced by Malkin ||l^. Theorems which deduce the stability of a nonlinear system from its 
linear first approximation usually (always?) assume the linear first approximation to be uniformly 
stable or uniformly asymptotically stable |^2| . The uniformity assumptions are not explicitly stated 
sometimes, for example in and In these cases, the A{t) in y(t) = A{t)y is either constant or 
periodic, which means that stability implies uniform stability and asymptotic stability implies uni- 
form asymptotic stability. Uniformity assumptions are natural in the theory of Lyapunov functions 
as well [^]. We will find them useful for bounding E{t). 

Definition 3. The solution x{t;xo) of x{t) = f{t,x) is uniformly stable if for every e > there 
exists a, 6 > such that ||x(r;xo) — a^rll < ^ for r > implies ||x(i;xo) — x{t;T, x'^)\\ < e for t > r. 



Definition 4. The solution x(t;xo) is uniformly asymptotically stable if it is uniformly stable and 
the choice of 5 in the previous definition can be made in such a way that \\x{t; xq) — x{t; t,x'^)\\ 
as r — > oo in a uniform way; i.e. given e' > there exists T^/ such that \\x{t;xQ) — x{t;T,x'^)\\ < e' 
for alH > r + T^/ and x'^ satisfying ||x(r; xq) — x'^\\ < 5. 



Theorem 5.5. // the zero solution of y{t) = A[t)y, y(0) = 0, is uniformly stable, its E{t) is 
linearly hounded; i.e., E{t) < Kt for some K > and < t < oo. 



Proof. Uniform stability of the zero solution is equivalent to boundedness of ||y(t)y "^{s) 
t > s > pTf. If \\Y{t)Y-'^{s)\\ <KfoTt>s>0, Corollary O implies E{t) < Kt. 



for 
□ 



Theorem 5.6. If the zero solution of y{t) = A[t)y, y{0) = 0, is uniformly asymptotically stable, 
its E{t) is hounded by a constant; i.e. E{t) < K for some K > {) and < t < oo. 

Proof. Uniform asymptotic stability of the zero solution is equivalent to ||y(t)y~"'^(s)|| < Me~'^''*~^^ 
for > 0, M > 0, and t > s > [^] |27|. Again, we can use Corollary |5.3| to complete the proof. □ 



Theorem 5^ implies that E{t) for the solution of x{t) = Ax, x(0) = xq, is linearly bounded 
if all the eigenvalues of A have negative or zero real parts, and the ones with zero real part are 
simple. If all the eigenvalues of A have strictly negative real parts then, in fact, E{t) is bounded 
by a constant by Theorem 5.6. The necessary stability properties of the zero solution of x{t) = Ax 



are verified in numerous places including |22|. 
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6 E{t) for Nonlinear Systems 



This section gives two approaches to the analysis of E{t) of nonlinear systems. Theorem 3^ proves 
that the E{t) of the solution of a nonlinear system and the E{t) of the zero solution of its first 
approximation are the same. So one approach is to look at the linearized problem. The other 
approach is to directly make stability assumptions about the solution of the nonlinear system 
(Theorems 6.5 and |6.6|). Both these approaches are illustrated in Section 7. 



Lemma 6.1. Assume that F(t,0) = 0, that F{t,y) has continuous first order partial derivatives 
with respect to y, and that F{t,y) is continuous with respect to t, t> 0. Then, 

F{t,y) = A{t)y + g{t,y), 

where A(t) = ^^^^|y_Q and g{t,y) = o(|[y||) as y ^ uniformly over compact intervals oft. 

Proof See [|. □ 

As noted in the introduction, the following theorem is implicit in the work of Henrici ||l^ |13|. 
It might be possible to generalize it to partial differential equations. 

Theorem 6.2. As in Lemma 6.1, let F{t,0) = 0, let F(t,y) have continuous first order partial 
derivatives in y and be continuous in t. Then the zero solution of 

yit)=F{t,y), y(0)=0, 

and the zero solution of 

y{t) = A{t)y, y(0) = 0, 
where A{t) = ^^gy^^ | have the same E(t). 

Proof. Let F{t,y) = A(t)y + g{t,y) as in Lemma 6.1. If needed, take h < hQ(t,r) for some r > 0. 
To reduce clutter in the proof, denote y{kh; h) by y^ and Y{kh) by Y^; Y{t) is the principal 
fundamental matrix of y{t) = A{t)y. Let n = \ t/h\ and = t — nh. 

Routine application of the variation of constants formula Q gives 

i-kh 

Vk = YkYk-iyk-i + / yky'^{s)g{s,y{s] h)) ds + ha{h)vk, 

J(k-l)h 



for k = 1,2, . . . ,n, and 



y{t;h)=Y{t)Y~'yn+ f Y{t)Y-\s)g{s,y{s;h)) ds. 

J nh 

Here, yo = and YkY-'^{s) = I + 0{h) for {k - l)h < s < kh. Further, \\y{s;h)\\ < te^^a{h) for 
< s < t by Lemma |3.1| , and = o{y) as y — > uniformly over s £ [0,t]. Therefore, we 

can write Yf:Y~^{s)g{s,y{s; h)) = Ca{h)r]{h)u{s) for < s < t, a constant C which can depend on 
t but not on h, and ||ti(s)|| < 1, where u{s) G W^; here ry(/i) — > as /i ^ 0. The expressions for 
and y{t; h) become, 

i-kh 

Vk = ykYk-iVk-i + Ca{h)r]{h) / u{s)ds + ha{h)vk, 

J(k-l)h 

y{t; h) = Y{t)Y-^yn + Ca{h)r]{h) I u{s)ds. 

J nh 
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Combining these equalities, 



^ = h Y{t)Y,7'vk + Chvih) Y{t)Y,-' r u{s)ds 

t 



+ Cr]{h) I u{s)ds. 

J nh 

Both the 2nd and 3rd terms vanish as, h ^ Q because of ri{h). Thus the nonlinear term g{t,y) has 
no effect on E(t). The proof can be formally completed using (|5.1[). □ 



Corollary 6.3. Let f{t, x) be continuous in t and have continuous first order partial derivatives 
with respect to x. The E(t) of the solution x{t;xo) of x{t) = f{t,x), 3;(0) = xq , and the E{t) of 
the zero solution of y(t) — A.(t)y, y(0) — 0, where A(t) — '^fe ^ \x—x{t-xo)' same. In fact, 



E{t) = snp [ ^^v{s)ds 
v(s) Jo dx{s) 



v{s) 

where the supremum is taken over continuous functions v{s) satisfying ||f (s)|| < 1. 

Proof. For a proof of the above formula for E{t), note that if Y{t) is the principal fundamental 
matrix of y{t) = A{t)y, then Y{t)Y~^{s) = for a proof see any basic book on ODEs. Plug 

this into Theorem to get the formula. □ 



For comparison with Theorem 3.2, we state Theorem 33 of Chapter IX of Sansone and Conti 



|22|. This result is typical of Lyapunov's theory of first approximation. 

Theorem 6.4. Assume the zero solution of y{t) = A{t)y, t >0, is uniformly asymptotically stable. 
If g{t,y) is continuous in t and y, and g{t,y) = o(||y||) uniformly over t > as y ^ 0, then the 
zero solution of y{t) = A(t)y + g(t, y) is also uniformly asymptotically stable. 

The assumption about g{t,y) being o(||y||) uniformly over the semi-infinite interval t > in 
Theorem 6.4 is quite stringent. A nonlinearity of the form g{t,y) = t[yf, . . . does not satisfy 



that assumption. But the theorem does not hold if the assumption is weakened to what is known 
about g{t, y) from Lemma 6.1; for a counterexample, see Bellman |l], p. 87]. A stringent assumption 
about (^(t, y) is not required in Theorem because how small an h we take to get convergence of 
e{t] h)/a{h) to E{t) can depend upon t. In definitions of stability, in contrast, we need one small 
perturbation which stays "close" till t = oo. 

The rest of this section is about bounding E{t) by making stability assumptions about the 
solution of a nonlinear system. For linear systems, this program was carried out in Section 4. We 
introduce a technique for error analysis of approximate solutions which uses Lyapunov functions. 
But first an example to point out some difficulties. 

Example 6.1. Let us first consider the zero solution of y{t) = y — e*y^, y{0) = 0, t > 0. Its 
E{t), by Theorems |6.2| and 5.1, is e* — 1. But we show that the zero solution is actually uniformly 



asymptotically stable. It is even exponentially asymptotically stable. 

Figure |^ is the portrait of trajectories of y(r) = y — e*y^. The portrait for ?/ < is a reflection 
about y = 0. Thus we can restrict ourselves to trajectories which are always in the upper half 
plane. Every point (t, y) with y > e~*/^ is on a trajectory that is pointed downwards. 
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Figure 3: The portrait of trajectories of y{t) = y — e y . All solutions tend towards the curve 
y{t) = e-*/2. 

We now verify uniform stability of the zero solution using the last observation. Given e > 0, 
choose so that e"**"/^ < e. By continuity properties, we can choose a 5 > so that if t < and 
y-r < 6, the trajectory through (r, y^) stays below e till t^. The maximum possible height (along y) 
of a trajectory beginning at (r, y,-), t > t,: and yr < S, is bounded by the larger of e""^/^ and \yr\- 
Since e""^/^ < e and 5 < e, uniform stability is verified. 

The verification of uniform asymptotic stability will be sketchy. We base it on the following 
facts: 

1. The solution of y{t) = y — e*y^, y(0) = yo tends to zero as i ^ oo for < yo ^ 1) 

2. Further, y{t; yo) < y{t; 1) for t > 0, < yo < 1, 

3. And, < y(t + r; r; yo) < y(t; 0; yo) for any yo > 0, r > 0, t > 0. 

The proof of item 1 involves a bit of elementary work which we do at the end. Item 2 is trivial. For 
item 3, think of y{t; 0, yo) as the solution of y{t) = y — e*y'^, y(0) = yo, and of y{t + r; r, yo) = z{t) 
as the solution of z{t) = z — e^'^'^z^, z{0) = yo, and use a differential inequality |jll|, p. 27]. 

Now let Te be such that \\y{t; yo)\\ < e for t > and yo = 1. Then \\y{t + r; r, yo)i| < e for any 
r > 0, t > Te and jyo| < 1. Thus, the zero solution is uniformly asymptotically stable. 

To prove item 1, it is enough to verify that the solution with the initial condition y(0) = 1 
satisfies y{t) < 2e~*/^ for t > 0. This is obvious from the portrait of trajectories: every trajectory 
that cuts the curve y = 2e~*/^ is in the downward direction; therefore any trajectory that starts 
below that curve has to stay below it for t > 0. In fact, y{t) < 2e~*/^ for the trajectory with 
y(0) = 1 implies that the zero solution is exponentially asymptotically stable, if we adopt Yoshizawa's 



definition of exponential asymptotic stability |27|. 



Here we have an example which is uniformly asymptotically stable, yet has an E(t) which 
increases exponentially with t. This is possible because the approximations y{t; h) can introduce 
errors at all points kh while the definitions of stability allow just one perturbation. 



The next two theorems are nonlinear analogues of Theorem 5^ and 5^. The proofs this time 
rely on the theory of Lyapunov functions. Let us note that Lyapunov functions V{t, y) are always 
assumed to be continuous in t and locally Lipshitz in x. Vp(t,y) is the rate of increase of V(t,y) 
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along a solution of y{t) = F{t, y) which goes through (t, y). More precisely, if y{t + r; t, y) is such 
a solution, then 

s y V{t + 6,yit + S-t,y))-V{t,y) 
ypit^y) =hmsup . 

If V^it, y) < aV{t, y) then V{t + 5, y{t + 6; y, t)) < e^'^Vit, y), and if V^(t, y)<0 then V{t + S, y{t + 
S\y,t)) < V{t,y); these two facts can be inferred from differential inequalities P [[Til]. 

We say that F{t,y) is uniformly Lipshitz in a neighbourhood of zero if \\F{t,yi) — F{t,y2)\\ < 
L\\yi ~ 2/211; for a constant L > and any yi with \\yi\\ < r, \\y2\\ < r where r > 0; L is the same 
constant for any t > 0. 

Theorem 6.5. Let F{t,y) be uniformly Lipshitz in y in a neighbourhood of 0. Assume that the 
zero solution of y{t) = F{t,y), y{0) = 0, is exponentially stable in the sense that 

\\y{t;to,yo)\\<Ke'<'-'^'^\\yo\\ 

for \\yo\\ < r, to > 0, t > to, c > 0, and K > 0. Then E[t) of the zero solution is bounded above by 
a constant. 

Proof. Stability assumptions in the theorem imply existence of Lyapunov function with following 
properties (Yoshizawa |2^, p. 97], Hale [1lO|] ): 

1- ||y|| < y{i-,y) ^ Cllylli where C > is a constant, 

2. \V{t,yi)-V{t,y2)\<L\\yi-y2l 

3. Vp{t,y) < —qcV{t,x) for some < q < 1. 

The domain of V{t, y) is t > and ||y|| < r for r > 0. 

Take h < h(){t,r). Because of item 3, V{t,y{t; h)) decreases at least by a factor e~'^'^^ along the 
approximate solution when t increases from kh to {k + on that interval of t the approximate 
solution follows the exact solution till the discontinuity at t = {k + l)h. The discontinuity can cause 
an increase in V{t,y) of at most L times its magnitude by item 2. Therefore, 

V{kh; y{kh; h)) < e'^'^^^Viik - l)h, y{{k - l)h; h)) + Lha{h), 

for A; = 1, 2, . . . , n and 

V{t,y{t-h)) < e-'"'''-V{nh,y{nh;h)). 
Combining these inequalities, we get 

/ 1 _ p—nqch . 

V{t,y{t-h))<e-'^^'^La{h){-^--^) 
< Ka{h). 

That K above can be a constant independent of h and n can be deduced from basic calculus using 
h < /iQ. Now, by item 1, \\y{t; h)\\ < Ca{h) implying a constant upper bound for E{t). □ 

The difficulty in proving Theorem 6.5 directly using the norm |[-|| is that when K > 1 the 
discretization error might actually be amplified by a factor greater than 1 over any given time step. 
Since we have to make the worst possible assumption at every time step, the final bound on E(t) 
obtained this way will actually be exponential in t when K > 1. The proof of Theorem 6.5 uses a 
carefully constructed Lyapunov function to get around this difficulty. 
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Theorem 6.6. Assume as in the previous theorem that F{t, y) is uniformly Lipshitz with respect 
to t in a neighbourhood of y = 0. If the zero solution of y{t) = F{t,y), y(0) = 0, t>0, is uniformly 
asymptotically stable, then E{t) < Kt for some constant K. 

Proof. Stability assumptions in this theorem imply the existence of a Lyapunov function with the 
following properties: (Hale ||lO| , Theorem 4.2, Chapter X], Yoshizawa ||27|| ) 

1. \\y\\<V{t,y), 

2. y(t,o) = 0, 

3. Vp{t,y) < 0, where Vp{t,y) is defined as in the previous proof, 

4. \V{t, yi) — V{t, 2/2)1 < K\\yi — 2/2!! for some constant K > 0. 

The domain of definition of V{t, y) is the same as in the previous proof. 
The proof is similar to that of Theorem 6.5, but this time 

V{kh, y{kh; h)) < V{{k - l)h, y{{k - l)h; h)) + Kha{h) 

for /c = 0, 1, . . . , n — 1 and 

V{t,y{t-h)) < V{nh,y{nh;h)). 

Combining these inequalities, we have V{t, y{t; h)) < Kta{h). As before, \\y{t; h)\\ < Kta{h), which 
this time implies that E{t) < Kt. □ 



Corollary 6.7. Let x{t;xQ) be the nonzero solution of the initial value problem x{t) = f{t,x), 
x{0) = xq. Assume f{t,x) is uniformly Lipshitz in x in a neighbourhood of x{t;xQ). Then uniform 
asymptotic stability of x{t;xo) implies that its E{t) is linearly bounded. 

Let us call attention to the necessity of making a uniform Lipshitz assumption about F{t, x) 



or f{t,x) in Theorem p.q or Corollary |6.7| . Example |6j, which is uniformly asymptotically stable, 
does not satisfy the uniform Lipshitz assumption and has an E{t) which increases exponentially. 
We do not know if Theorem |6.6| is still true if the assumption of uniform asymptotic stability is 
weakened to just uniform stability. If such a theorem were true, its wider applicability might be of 
use. Table ffl summarizes all of Sections 5 and 6 except Theorem S.2. 



7 Three Applications to Dynamical Systems 

We give three examples to illustrate the applicability of our methods for bounding the accumulation 
of global error. 



7.1 Hyperbolic Sinks of Dynamical Systems 



Let p be a fixed point of a dynamical system x{t) = f{x); i.e. let f{p) = 0. Then p is a 
hyperbolic sink, if all the eigenvalues of §j|^_ have strictly negative real parts. The following 
theorem can be derived using Chapter 6 of |25]. We give the theorem here because our method of 
proof is different. 



Theorem 7.1. Letx{t;xo) be a trajectory of the dynamical system x{t) = f{x), f G C^{R'^), which 
falls into a hyperbolic sink p as t ^ 00. Then its E[t) is bounded above by a constant. 
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Linear Problems 


Stability 


E{t) can increase exponentially 
even if ^(t) is bounded 


Asymptotic stability 


E{t) can increase exponentially 
even if is bounded; Exam- 
ple 


Uniform stability 


E{t) must be linearly bounded; 
Theorem ^.5| 


Uniform asymptotic stability 


E(t) must be bounded by a con- 
stant; Theorem |5.6| 


Nonlinear problems 


Uniform stability 


E{t) can increase exponentially 


Uniform stability with uniform 
Lipshitz assumption 


Not known if E{t) must be linearly 
bounded 


Uniform asymptotic stability 


E{t) can increase exponentially; 
Example |6.1| 


Uniform asymptotic stability with 
uniform Lipshitz assumption 


E{t) must be linearly bounded; 
Theorem |6.6| 


Exponential stability as in Theo- 
rem 6.5 


Not known if E{t) must be linearly 
bounded 


Exponential stability as in Theo- 
rem 6.5 with uniform Lipshitz as- 
sumption 


E{t) must be bounded by a con- 
stant; Theorem 6.5 



Table 1: Summary of part of Sections 5 and 6. The second column is the stability assumption 
about the solution; the last column says what is known about the conditioning function E{t) 
corresponding to such a solution. 
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Proof. Without loss of generality, take p = 0. By [20, p. 150], there is a neighbourhood Uq of 
such that xq G Uq implies 



||x(t;xo)|| < ce 

for constants a > 0, c > 0, and for t >0. 

Using continuity properties of solutions of differential equations (or openness of the basin of 
attraction), we can assume an open neighbourhood U of {x{t;xQ)\t > 0} in R'^ and a compact set 
K containing U such that K and consequently U are both contained in the basin of attraction of 
p = 0. Since ah trajectories beginning in K enter Uq in a finite amount of time, which depends 
only on iC, we can assume 

||x(t;xo)|| < Ce-"^ 

for constants a > 0, C > 0, for t > 0, and for any xq € K. We can also take K and U to be 
invariant under the flow. 

For any trajectory x{t; xq), i > 0, with xq € U , there exists r > such that if ||x(t; xq) — xi || < r, 
t > then xi S K. For such an xq, if < r, 

||x(t;xo) - x(t;r,x(r;xo) + S)\\ < 2Ce~"-^ 

for t > T > 0. If Xq fiJ but is in the basin of attraction of p, the trajectory x(t; xq) enters U in 
a finite amount of time. So we can get the same kind of bound as above by adjusting C and a if 
necessary. 

To apply Theorem |5.5| and deduce boundednes of E{t), we need only verify uniformity of the 
Lipshitz condition on /(x) for x £ K. This is trivial since K is compact and /(x) is C^. 

□ 

7.2 Hyperbolic, Attracting Cycles of C^+'' Dynamical Systems 

Let x(t), t > 0, be a periodic orbit of the dynamical system x(t) = /(x) in R'^. Let T > be 
its period so that x{t + T) = x(t). Denote the set of points on this orbit by 7. 

The characteristic multipliers of the cycle 7 can be defined in two ways. One is to pick a point 
p € 7, take a cross-section E at p, define a Poincare map for S, and then define the characteristic 
multipliers as the {d — 1) eigenvalues of the linearization of the Poincare map at p. The other way is 
to consider the linear first approximation y{t) = A{t)y on the cycle 7. Obviously, A{t + T) = A{t) 
for t > 0. The Floquet numbers of this linear system can also be used to define characteristic 
multipliers. For a lucid account of these matters, see Robinson |20|] . 

The cycle 7 is hyperbolic and attracting if all its characteristic multipliers are strictly less than 
1 in magnitude. 

Theorem 7.2. Let x{t; xq), t > 0, be an orbit of a C^"*"*^ dynamical system x{t) = /(x) in R"^ which 
falls into a hyperbolic, attracting cycle 7 as t ^ 00. Then its E{t) is linearly bounded from above. 

Let us first prove the following lemma. If any one solution of a linear system is uniformly stable, 
so is every other solution. So we might speak of the linear system itself as being uniformly stable. 

Lemma 7.3. Assume xq £ ^ so that x(t;xo) is a periodic orbit. Let its linear first approximation 
be y{t) = A{t)y, t >0. If j is hyperbolic and attracting, y{t) = A{t)y is uniformly stable, and the 
E{t) associated with x(t;xo) is linearly bounded. 



21 



Proof. Uniform stability of y{t) = A{t)y is an easy consequence of the characteristic multipliers 
of 7 being strictly less than 1. See Chapter IX of [22|. The linear bound on E{t) is implied by 

□ 



Theorem |5.5| and Corollary 6.7. 

Lemma is contained in a different form in the work of Cano and Sanz-Serna |^]. But let 
us note that Theorem 7^ goes beyond Lemma 7^3 in a significant way. In practice, it is highly 
unlikely that xq itself is on the cycle 7. But it is often easy to find xq so that x{t;xQ) falls into a 
cycle 7. 

The following lemma is known as the Dini-Hukuhara-Caligo theorem. It is Corollary 1 of 
Chapter IX of [^]. Its proof, which we omit, is short and simple, and illustrative of an important 
technique in stability theory. 

Lemma 7.4. Assume the linear system y{t) = A{t)y, t >0, is uniformly stable. Assume also that 
B{t), t > 0, is continuous with J^\\B{t)\\dt < 00. Then the linear system y{t) = {A(t) + B(t))y is 
also uniformly stable. 



Proof of Theorem ''Li. By Hartman |11, p. 254], there exists a point Xq G 7 such that 

||a;(t;a;o) - x(t;xo)|| < ce""*. 



(7.1) 



for constants a > 0, c > 0, and for t > 0. This is called convergence in phase [20|. 

Let y(t) = where ^(t) = be the first approximation along x(t;xo). By 

Lemma [7.3| , this linear system is uniformly stable. 

Let y(t) = {A{t) + B{t))y, where A{t) + B{t) = ^\^^x{t-xo) ^'"^^ approximation along 

x{t;xo). The estimate (0) for convergence in phase implies 



< cie 



-ai t 



for constant ai > and ci > 0. This is because both x{t;xQ) and x{t;xQ) stay within a compact 
region of R'^, and f{x) is C^^*^ (this where we need the e in C^"*"*^). By Lemma 7.4, y{t) = 
{A{t) + B{t))y is also uniformly stable. 

Since A{t) + B{t) gives the linearization of x{t; xq), the proof is easy to complete using Theorem 
5.5 and Corollary |6.7|. □ 



7.3 Normally Contracting Manifolds with Quasiperiodic Flows 

Let us introduce the notation 0j for the flow induced on R'^ by x{t) = f{x). With this notation 
x{t;xo) = 4>tXo. Let F be a compact manifold which is invariant under this flow. We consider 
the situation when the flow on V is differentiably conjugate to quasiperiodic flow on a torus, and V 
is normally hyperbolic and contracting, or briefly, normally contracting. We now explain the two 
italicized concepts in this paragraph. 

A torus T" is n copies of the circle . If the angle on the ith circle is parameterized by ^j, a 
quasiperiodic flow on T" is of the form 6i{t) = {Oi{0) + ait) mod 27r. In fact, the flow is periodic 
if the Ui are all mutually commensurable. Denote this flow by ij^t- 

When we say that the flow (j)t is differentiably conjugate to quasiperiodic flow on a torus, we 
mean that there exists a homeomorphism h :V ^ T" such that h{(j)tx) = ■0i(^x) for x & V and 
t > 0. 
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To define normal contractivity |2C], associate a direct sum decomposition Nx of 
with every x in V. In this sphtting is the tangent space of V at x, and Nx, the normal space, 
varies continuously with x. If the Nx can be chosen so that 

||niv,|^|A^x|| <ce-'^*, 

where y = (j)tx, the matrix inside the norm is the restriction of the derivative ^ to act from Nx to 
Ny, c > 0, and yu > 0, then V is normally contracting. Usually, the definition of normal contractivity 
comes with an other assumption which says contraction in the normal direction dominates any 
contraction on the manifold V. But since we have assumed that the flow on V is differentiably 
conjugate to quasiperiodic flow on a torus, this other assumption can be dropped. 

Obviously, the tangent spaces Tx are invariant under the derivative map ^^p- It is actually 
possible to choose Nx so that they too are invariant under the derivative map |14] [pH]. We take 



this to be the case. So the derivative map maps Tx to Tfp^x and Nx to N^^x- 

Theorem 7.5. Let x{t; xq) be a trajectory of the C^^*^ flow x{t) = f{x) which falls into a normally 
contracting and invariant manifold V. Assume that the flow on V is differentiably conjugate to 
quasiperiodic flow on a torus. Then the E(t) of x{t;xo) is linearly bounded. 



The above theorem generalizes Theorem |7.2| . Its proof is exactly analogous. We begin with a 
lemma about trajectories that begin on V. 

Lemma 7.6. Let xq (^V so that x{t;xo) stays on V for t>0. Make the same assumptions about 
V as in Theorem [7. 4 - Then the linearization y{t) = A{t)y along x(t;xo) is uniformly stable, and 
the E{t) of x{t;xo) is linearly bounded. 

Proof. The principal fundamental matrix of the linearization in the lemma is given by Y{t) = 

which is the derivative map. Therefore, it is enough if we show that || || is bounded by a constant 
for any x GV and t > 0. 

We already know that the maps induced by the derivative map between tangent spaces and 
between normal spaces are bounded in norm because of differentiable conjugacy to flow on a torus 
and normal contractivity, respectively. Since the tangent spaces and normal spaces are both invari- 
ant under the derivative map, it is enough if we show that the angle between Tx and Nx (in the 
sense of the CS decomposition) is bounded away from 0. That this angle is bounded away from 
is implied by the compactness of y. □ 



The proof of Theorem 7.5 can be completed exactly like the proof of Theorem 7.2 using the 



following result about convergence in phase. 

Theorem 7.7. As in Theorem \7.^ let x{t;xo) be a trajectory of the flow x{t) = f{x) which 
falls into a normally contracting and invariant manifold V, and let the flow on V be differentiably 
conjugate to quasiperiodic flow on a torus. Then there exists x'q £V such that 

\\x{t;xo) - x{t;xQ)\\ < ce~'**, 

for positive constants c and a, and t > 0. 

Proof. This theorem can be deduced from Theorem 4.1 and the remark following its proof in Hirsch, 



Pugh and Shub |14|. See in particular part (a) of that theorem about stable manifolds and part 



(g) about conjugacy to linearized flows. □ 
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8 Three Remarks 



This last section is a collection of parenthetical remarks about matters which are related to E{t) 
and which we have not investigated in detail. 

(i) Multistep methods and variable time stepping. The model for discretization errors in Section 2 
is adapted to single step methods with constant step sizes. For linear multistep methods with 
constant step sizes, we believe the accumulation of global error can be worse but not better 
(after excluding some trivial cases) than indicated by E{t). 

For variable time stepping, it is sometimes but not always true that the ratio of the largest 
to the smallest time step is bounded by a constant p3| [24|. In this case at least, the effect 



of variable time stepping is to improve global errors by a constant factor over the indication 
given by E{t). 

(ii) One sided Lipshitz conditions. Let us briefly summarize the approach to global error analysis 
using one sided Lipshitz conditions for the linear system y{t) = A{t)y. Since = 
y^{t)y{t), we have 

'^^^ = y^{t){A^{t)+A{t))y{t). 

If \{t) is the maximum eigenvalue of ^ {t) + A{t), then 

d\W)\? 



dt 



<m\\ym 



Thus upper bounds for ||y(t)||, and hence for ||y(t)|| where Y{t) is the principal fundamental 
matrix of the linear system, can be written down in terms of \{t). These can be plugged 



into Theorem 5^ to get bounds on E{t) and hence the accumulation of global error. For a 
detailed account, see 

In our view, one sided Lipshitz conditions are basically a way to get a handle on stability by 
looking at the evolution of the norm of y{t). This is a far less general approach to stability 
than the two methods of Lyapunov we have used (this deficiency of norms is something 
Lyapunov must have realized). It is also of far lesser applicability; we do not see a way to 
derive any of the results in Section 7 by using one sided Lipshitz conditions. 

(iii) Numerical symplectic integrators and Hamiltonian problems. 

Special methods which preserve properties like symplecticity or energy conservation or volume 
conservation in phase space show milder increase of global errors with t than general methods 
in numerically integrating Hamiltonian problems. Calvo and Sanz-Serna |^] showed using 
careful numerical experiments that the accumulation of global error with t was linear for a 
symplectic method but quadratic for a Runge-Kutta method when the methods were applied 
to Kepler's problem. A full explanation came after ^ and in the paper [Q] by Cano and 
Sanz-Serna. There are numerical experiments in Quispell and Dyt |19] which are unexplained. 



The conditioning function E{t) is not suitable for studying the accumulation of global error in 
these special methods. The class of local discretization errors allowed by the model in Section 
2 is too general. The discretization errors made by these special methods are restricted 
in some way; for example, the discretization errors of an energy preserving method cannot 
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take the approximation out of the constant energy manifold. But is it possible to define 
conditioning functions which are adapted to these special methods? We think it might be 
possible. The concepts of stability needed to bound such conditioning fuctions will also, no 
doubt, be different. 
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